Improve handling of large values by inverse hyperbolic and trigonometric functions#2928
Improve handling of large values by inverse hyperbolic and trigonometric functions#2928antonwolfy wants to merge 12 commits into
Conversation
…x_math::cacos() function
|
View rendered docs @ https://intelpython.github.io/dpnp/pull/2928/index.html |
|
Array API standard conformance tests for dpnp=0.21.0dev0=py313h509198e_54 ran successfully. |
should this be filed as an issue with DPC++ compiler? Since most implementations should have logic to prevent overflow/infs in such cases |
Yes, we probably should, but I wonder that the extension still in the experimental state. I also tried to implement the kernel through SYCL's C complex math functions since marked as supported extension: // Forward declare SYCL's C complex math functions with C linkage
extern "C" {
extern __DPCPP_SYCL_EXTERNAL_LIBC float __complex__ cacosf(float __complex__ z);
extern __DPCPP_SYCL_EXTERNAL_LIBC double __complex__ cacos(double __complex__ z);
}
...
if constexpr (is_complex<argT>::value) {
using value_type = typename argT::value_type;
if constexpr (std::is_same_v<value_type, float>) {
auto result = ::cacosf(*reinterpret_cast<const float __complex__*>(&in));
return *reinterpret_cast<const std::complex<float>*>(&result);
}
else {
auto result = ::cacos(*reinterpret_cast<const double __complex__*>(&in));
return *reinterpret_cast<const std::complex<double>*>(&result);
}
}
...but faced exactly the same #2924 issue, since the SYCL does provide double __complex__ w = clog(z + csqrt(__sqr(z) - 1.0));Btw I'd propose to merge that PR, in background I'll submit the issue towards the compiler. |
The
dpnp.tensor.acoshfunction was returning incorrect results (infinity) for complex numbers with large negative real parts (see #2924).The issue was caused by the threshold used
1/epsilonto decide when to switch from acos in SYCL's experimental complex extension (exprm_ns::acos()) to a log-based formula.The SYCL
exprm_ns::acos()returns incorrect results for large negative real parts, due to used formulaz + sqrt(z² - 1), which leads tolog(0) = -infinity.The PR proposes to changed the threshold from
1/epsilontosqrt(1/epsilon)/2for all inverse trig/hyperbolic functions. This is the precision loss point wheresqrt(z² ± 1)calculations become unstable.Also the PR reduces code duplication by creating shared
casinh(),casin(),catanh(),catan()functions in newcomplex_math.hppheader.